Mean Field Fluid Behavior of the Gaussian Core Model 
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We show that the Gaussian core model of particles interacting via a penetrable repulsive Gaussian 
potential, first considered by Stillinger (J. Chem. Phys. 65, 3968 (1976)), behaves like a weakly 
correlated "mean field fluid" over a surprisingly wide density and temperature range. In the bulk 
the structure of the fluid phase is accurately described by the random phase approximation for the 
direct correlation function, and by the more sophisticated HNC integral equation. The resulting 
pressure deviates very little from a simple, mean-field like, quadratic form in the density, while the 
low density virial expansion turns out to have an extremely small radius of convergence. Density 
profiles near a hard wall are also very accurately described by the corresponding mean-field free- 
energy functional. The binary version of the model exhibits a spinodal instability against de-mixing 
at high densities. Possible implications for semi-dilute polymer solutions are discussed. 
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I. INTRO 

Interactions between atoms or molecules in simple flu- 
ids invariably contain a short-range repulsive component 
or hard core, such that the local molecular structure is 
dominated by excluded volume effects. This observation 
explains the success of simple models involving hard con- 
vex bodies in explaining the structure and phase transi- 
tions in simple atomic or molecular fluids 0. For ex- 
ample, the hard sphere model has been instrumental in 
understanding freezing of simple fluids . The same ex- 
tends success to somewhat more complex fluids like liquid 
crystals, where hard ellipsoids or spherocylinders have 
been widely used to investigate the isotropic-to-nematic 
transition and other mesophases [||. However the situa- 
tion is generally not as simple in complex fluids, where 
interactions between mesoscopic particles are often of ef- 
fective and of entropic origin. While excluded volume 
effects still dominate the interaction between compact 
colloidal particles, the effective forces between "soft" or 
fractal objects of fluctuating shape, like polymer coils 
or membranes, cannot be modeled by hard cores. Poly- 
mers in a good solvent form highly penetrable coils and 
it is by now well established that the effective interaction 
between the centers of mass of two polymer coils, duly 
averaged over internal conformations, is finite for all dis- 
tances, and decays rapidly beyond the radius of gyration 
of the coils . For two isolated non-intersecting poly- 
mer chains, the effective pair potential at zero separation 
of the centers of mass v(r=0), is of the order of 2ksT 
for sufficiently long chains |6|7|], and is reasonably well 
represented by a Gaussian whose width is of order the 
polymer radius of gyration Rq, as shown in Fig. [l]. 

We have recently shown that the general shape of the 
effective pair potential remains roughly the same in dilute 
and semi-dilute solutions of self-avoiding random walk 
(SAW) polymers, and does not vary strongly with poly- 
mer concentration (cf. Fig. [I]) ||. The effective pair po- 
tential model has been shown to accurately reproduce the 



structure and thermodynamics calculated from Monte 
Carlo (MC) simulations of solutions of SAW polymers 
over a wide range of concentrations 0]. 
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FIG. 1. Polymer center of mass potentials /3v(r) 
from simulations of L = 500 monomer SAW chains 
Q| are compared to a best-fit Gaussian ([!]), deter- 
mined by fitting (3v(0) to fix /3e, and f3v(0) to fix 
R. The potential for two isolated coils (p — * 0) 
is well approximated by a Gaussian potential with 
f3e — 1.87, R = 1.13-Rg- The potential in the semi-dilute 
regime (p ~ 4x 3/(47r_Rg)) is approximated by a Gaussian 
potential with (3e = 2.16, R = l.45R G - 



Neglecting in the first instance the state dependence 
of the effective potential, it seems hence worthwhile to 
examine the equilibrium properties of a fluid of "soft" 
particles interacting via a pair potential approximated 
by a simple Gaussian form: 



v(r) = eexp 



R 2 



(1) 



1 



where e is the energy scale, and R determines the width. 
The Fourier transform (FT) is 

v(k) = 7r 3 / 2 i? 3 e exp(-^) . (2) 

Such a "Gaussian core model" (GCM) was in fact intro- 
duced some time ago by Stillinger M, who focussed on 
the low temperature regime e* = e/ksT >> 1, where 
the model exhibits hard-sphere like behavior, and a re- 
entrant fluid-solid-fluid phase diagram under compres- 
sion below a threshold temperature. This work was fur- 
ther expanded by Likos et al. |n|, who showed that the 
model remains fluid at all densities when e* < 100. They 
also demonstrated that for this model, the familiar HNC 
closure for the pair distribution function g(r) becomes 
exact in the high density limit, and that the random 
phase approximation (RPA) is remarkably accurate at 
high densities. 

In this paper we concentrate on the fluid phase of the 
GCM (e* < 100), with a particular emphasis on the 
regime relevant for polymer solutions (e* ~ 2) S, for 
which the dilute regime corresponds to reduced densities 
p* = pR 3 < 3/(47r) w 0.239, and the semi-dilute regime 
corresponds to p* > 3/(4n) @ (here p = N/V is the 
number of Gaussian core particles per unit volume). We 
shall successively consider the homogeneous fluid phase, 
the inhomogeneous fluid phase in the vicinity of a hard 
wall, and the possibility of de-mixing of binary Gaussian 
core systems. 



II. THE HOMOGENEOUS FLUID PHASE 

A. The Thermodynamic stability of the GCM fluid 

We consider a system of N particles interacting via a 
Gaussian pair potential (|I|), in a volume V. In the ab- 
sence of an infinitely repulsive core the first question is 
that of thermodynamic stability against collapse, i.e. the 
existence of a well defined thermodynamic limit. Accord- 
ing to definition 3.2.1. in Ruellc's classic book [j^J, the 
total interaction energy Vn, which can be built up of 
pair and higher order potentials, is stable if there exists 
a B > such that 

F/v(ri, r N ) > —NB (3) 

for all N > and all {n} in the phase-space R N . Sta- 
bility implies convergence of the grand partition function 
and a well defined thermodynamic limit. Specializing to 
pair potentials i>2 , the total potential energy of the sys- 
tem, for any configuration of N particles {ri} £ R N , can 
be written as : 

4 2) ( ri ,....,r N )= J2 Vaflri-rjl) (4) 

l<i<j<N 



For purely repulsive pair potentials, such as the GCM 
with e* > 0, satisfies the condition (||), so that a well 
defined thermodynamic limit exists. However, if V2(r) is 
not strictly positive, this may no longer be true. In Ap- 
pendix A two examples are discussed, involving a finite 
core and (small) attractive tail, which do not lead to a 
proper thermodynamic limit. 



B. The Structure of the GCM fluid 

To determine the pair structure of the GCM fluid, we 
have used the HNC closure which becomes exact in the 
high density limit; this closure relates the direct correla- 
tion function c(r) to the pair potential v(r) and the pair 
correlation function h(r) = g(r) — 1, according to: 

c(r) = -0v{r) + h(r) - In [1 + h(r)} , (5) 

where = 1/ksT. This closure must be combined with 
the Ornstein Zernicke (OZ) relation between c(r) and 
h(r) fl3| ] to yield a non-linear integral equation, which 
must be solved numerically. Examples for e* = 3 at three 
reduced densities p* are shown in Fig |^, and compared 
to the results of MC simulations. 
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FIG. 2. Comparison of MC simulations and solu- 
tions of the HNC integral equation in a regime relevant 
for polymer solutions g|, v(r) = 2exp [-(r/R) 2 ]. The 
lines are HNC calculations, and the symbols represent 
MC simulations for different reduced densities p* . 

The key feature is that the "soft" correlation hole is 
gradually reduced as p* increases, a behavior typical of 
finite core potentials, which leads to overlap and ideal- 
gas like behavior of g(r) in the high density limit. Note 
that the HNC results are indistinguishable from the MC 
data, so that for e* ~ 2 the HNC correlation function will 
henceforth be considered as providing an "exact" refer- 
ence to gauge simpler theories. The simplest is the RPA 
p3|,p^ |, which may be formally derived from the HNC 
closure (^) by linearizing the logarithm, leading to: 



2 



c(r) = —j3v{r). 



(6) 



0.5 - 



Since Fig. g clearly shows that the amplitude of h(r) is 
rather small at high densities, we may expect the RPA 
closure to become more accurate as the density increases. 
For the GCM, Eq. (@) and Eq. (|) imply the following 
FT of c(r): 



c(k) = -e*7r 3/2 i? 3 exp 



-k 2 R 2 



(7) 



and the OZ relation immediately yields the following 
RPA structure factor: 



S(k) = 1 + ph(k) = 



1 



1 - pc{k) 



1 + acxp [-fc 2 i? 2 /4] ' 



(8) 



where we have introduced the dimensionless coupling pa- 
rameter: 



3/2 



(9) 



HNC results for c(r) and S(k) at several densities are 
compared to the RPA predictions in Figs. [| and ||. 

Since h(r) > ln[l + h(r)], the HNC direct correlation 
functions are bounded below by the RPA form (^|). Fig. 

also shows that the HNC c(r) appears to be bounded 
above by the low density approximation: 
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c{r) = f(r) = exp [— f}v(r)\ — 1 



(10) 



which corresponds to the lowest order term in the expan- 
sion of c(r) in powers of p p3| ; f(r) is the usual Mayer 
/-function. Figs. || and ^ also illustrate the point that the 
simple RPA becomes very accurate at high densities, so 
that it is worthwhile to inquire about a correction to Eq. 
(||). Expanding the logarithm on the r.h.s. of Eq. (||) to 
second order in h, one arrives at the following expression 
for c(r): 



c(r) = -j3v{r) + -h(rf 



(11) 



Solution of the closure (|ll|) and the corresponding OZ 
relation requires an iterative procedure, as for the full 
HNC closure. Further simplification amounts to replac- 
ing h(r) in Eq. ([ll]) by its RPA form derived from Eq. 
(§) by FT; we refer to this non-iterative approximation 
as RPA2 pi: 



c(r) = -0v(r) + hi R p A {rf 



(12) 



From Fig. |^ it is clear that crpa2 (r) is indistinguishable 
from the HNC results except at low densities (p* < 0.2). 
The limitations of RPA theory at low densities become 
apparent by considering the resulting behavior of g(r) at 
short distance. 
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FIG. 3. HNC, RPA, and RPA2 forms of the di- 
rect correlation function c(r) for v(r) — 2exp[— (r/R) 2 ]. 
From top to bottom the densities are p* = 0.1,0.5, 
and 1 respectively. Note that the HNC c(r) is 
bounded by crpa(j") = —/3v(r) from below and 
f(r) = exp[— f3v(r)] — 1 from above. Inset: Ratio's 
crpa(0)/c H nc(0) (solid line) and chp^^/chat^O) 
(long-dashed line) v.s. density p*. For p* > 0.05, 
crpa2(0) is always within 2% of chnc(0) 




FIG. 4. RPA and HNC forms of the structure factor 
S(k) for v(r) = 2 exp[-(r/i?) 2 ]. From top to bottom the 
densities are p* = 0.01, 0.1, 0.5, 5 respectively. 

The zero separation value is easily derived from the 
-> limit of the FT of Eq. (§), with the result : 



9rpa(0) = H Li 3/2 (-a), 

a 



where the nth polylogarithm is defined by: 

OO 

Li n (x) = J2 T^- 



(13) 



k=l 



(14) 
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for \x\ < 1. If e* < 1, gRPA(0) is positive for all den- 
sities p*. However, when e* > 1 there is always a re- 
duced density pg below which gppA{0) < 0, which is 
unphysical. For example, if e* = 2, gppA(0) < 0, for 
p* < Pq = 0.3617. However, even for p* < p* , the struc- 
ture factor S(k) is still reasonably well described by the 
RPA because the deficiencies of g(r) at small r do not 
strongly affect S(k). This is also illustrated in the inset 
of Fig. ||, where crpa2{0) is seen approximate the quasi- 
exact HNC result to within 2% for densities p* > 0.05, 
which is a significantly lower bound than = 0.3671. 



C. Thermodynamics of the GCM fluid 

Turning now to thermodynamic properties, the equa- 
tion of state can be calculated via either of two routes 
p3[ : from the compressibility equation: 



0P = 



p dpP(p') 
dp' 



dp' = [ P [l-p'c(k = 0;p')]dp' (15) 
Jo 



where P denotes the pressure. For soft-core potential 
systems the virial equation leads to: 



1 



/3P = p+y«(fc = 0) 



2tt 



r 6 „ w h(r)dr. 



dr 



(16) 
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FIG. 5. Compressibility factors (Z — flP/p) from 
RPA and HNC for a Gaussian potential with e* = 2, 
compared to MC simulations and to two and three term 
virial expansions. Z PPA , Zf[ NC , and Z V HNC , are indis- 
tinguishable on this scale. Inset: The analytic ratio 
Zrpa/Z v rpa = (l-e*H(a)/(l + a/2))- 1 gives a good ap- 
proximation to Z PPA /Zff NC and goes to zero as p — > oo. 
The ratio Z\ PA IZ V HNC demonstrates that Z PPA approx- 
imates the true e.o.s. to better than 1% accuracy over the 
entire density range for e* = 2. 



If the correlation functions h(r) and c(r) were known 
exactly, the two routes would lead to identical equations 
of state. Approximate theories are not, in general, ther- 
modynamically consistent. However, as shown in Figs. 
HI?], the HNC closure yields practically identical values 
of the pressure over the whole range of densities, even 
for a repulsion as high as e* = 90 (remember that for 
e* > 100 re-entrant crystallization sets in ||,[l(|). More- 
over, the two HNC estimates of the pressure agree closely 
with the results of MC simulations. These results con- 
firm the conjecture that HNC and RPA become exact for 
the GCM in the high density limit, but also show that 
the HNC works well for low densities in the fluid regime 
we consider (e* < 100). 




FIG. 6. Compressibility factors from RPA and HNC 
for a Gaussian potential with e* = 10, compared to MC 
simulations and to two and three term virial expansions. 
Z PPA , Zx NC , and Z V HNC are very close over much of 
the density range. Inset: Compressibility factors at low 
density, symbols are the same as in the main figure. Note 
that Z PPA shows unphysical behavior for very small p* , 
which can be understood from the effective virial expan- 
sion discussed in Appendix B. 



Turning now to the much simpler RPA, it is easily ver- 
ified from Equations (||) and ([l5]) that the dimensionless 
equation of state (e.o.s.), Z = (3P/p, reduces, within the 
compressibility route, to the simple expression: 



Irpa = 1 + 2 P ^ k = °) 



which for the GCM leads to: 



1 



Z RPA = 1 + 



(17) 



(18) 



This in turn leads to an excess free energy per particle: 
E ir = 2^(fc = 0)= (19) 



4 



identical to that obtained from a van der Waals like mean 
field theory (MF) so that Z C RPA = Zmf', it also implies 
that the excess chemical potential is linear in density. 



where: 
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FIG. 7. Compressibility factors from RPA and HNC 
for a Gaussian potential with e* = 90, compared to 
MC simulations and to two and three term virial expan- 
sions. Again, the HNC virial approximation is nearly 
exact across the whole density range, and the e.o.s. is 
that of a mean-field fluid at all but the lowest densities. 
The low-density limit is further discussed in Appendix B 
and illustrated in Fig. [L5|. 

Remembering that the quasi-exact HNC direct corre- 
lation function is bounded below by the RPA form (||), 
it is immediately clear from the compressibility equation 
(|l5|) that one may expect: 



Z < Zrpa — Zmf 



(20) 



This conjecture is supported by Eq. (fl6[), which shows 
that the exact equation of state is given by: 



z-z c 2 V 

^ — RPA TP 

JO 



r 6 „ ' h(r)dr, 



Or 



which, for the Gaussian core model reduces to: 
4a 



Z 



1 

l + -a 



3v^ Jo 



x e x h(x)dx 



(21) 



(22) 



where the dimensionless spacing x = r/R was intro- 
duced. The conjecture ( p0| ) is thus true, provided the 
integral on the r.h.s. of Eq. (|2^) is negative. This is very 
likely for sufficiently high temperatures since the HNC 
results plotted in Fig. || show that h(x) is mostly nega- 
tive. 



The integral in Eq. (22) can be calculated analytically 
within the RPA, leading to the following result for the 
RPA virial equation of state: 



Zrpa = 1 + ~<* - e*N(a) 



(23) 



2a 



Lis (— a) — Lis(—a) 



(24) 



K(a) is zero for a = 0, has a maximum of 0.0908 at 
a = 7.8, and goes to zero for a — > oo, which implies that 
for any e*, the RPA becomes thermodynamically consis- 
tent in the high density limit. 

Interestingly, within the RPA2, the compressibility 
e.o.s. may also be solved analytically and yields exactly 
the same result (j23|) as the virial e.o.s. in the RPA (i.e. 



J RPA2 



Zrpa) 



suggesting that the latter is more ac- 
curate than the MF or RPA compressibility equation of 
state (17). In fact, as shown in Figures the RPA 
virial equation of state is virtually indistinguishable from 
the practically self-consistent HNC results and the MC 
simulations, except at very low reduced densities p* . Fig. 
[?] demonstrates that the mean-field approximation ( |l~7| ) is 
still surprisingly good, even for an interaction as large as 
e* = 90, just below the value where freezing sets in. This 
implies that the hard-sphere limit, envisioned by Still- 
inger 0], is still not reached at such a strong interaction, 
and that the Gaussian core-model behaves very much like 
a "mean field fluid" (MFF) over a wide temperature and 
density range. 

Perhaps the most striking result is the persistence of 
the linear slope of the equation of state to such low den- 
sities |16|. The slope differs, however, from that deter- 
mined by the (smaller) second virial coefficient. This is 
further discussed in Appendix B, where it is shown that 
the standard virial expansion of the equation of state M] 
has a very small radius of convergence for the GCM, and 
is of limited use for this model, contrarily to the case of 
the hard sphere fluid p7|. 



III. THE GCM NEAR A WALL 

In view of the success of the HNC and RPA theories 
for the GCM in the homogeneous bulk fluid phase, it is of 
interest to ascertain the validity of these approximations 
under inhomogeneous conditions, e.g. in the presence of 
an external potential <fi(r) acting on the particles. In 
this section we shall consider more specifically the den- 
sity profiles of GCM particles near a hard wall, using the 
formalism of density functional theory (DFT). In an ex- 
ternal potential the density of particles will change from 
a constant bulk value, pb, to a spatially varying local den- 
sity p(r). The grand potential of the inhomogeneous fluid 
in equilibrium with a bulk reservoir fixing the chemical 
potential p, may be cast in the generic form Jl8[ 

/3fUp(r)] = PT m [p(v)] - J dv (Pp ~ 0<f>(r)) p(v) (25) 

where the intrinsic free energy functional T m naturally 
splits into ideal and excess parts, T %d and T ex . The lat- 
ter is an unknown functional of the local density p(r). 
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When the inhomogeneity is not too strong, the excess 
part T ex may be expanded in a functional Taylor series 
in the deviation of the local density p(r) from the bulk 
density p b . If the expansion is truncated after second 
order the HNC functional results: 



pn v [p(r)} = [3n[p b ] + J dr'/30(r')p(r') 
+ / dr'{p{v')\nW)/p b ]-p{v') + Pb } 



--J drdr'(p(r) Pb )c b 2 \\r r'\)(p(r') - Pb ). (26) 

This functional is to be minimized with respect to p(j), 
and the resulting Euler-Lagrange equation reads: 



p(r) = p b exp 



-/¥(r) 



Pb 



dr'c 



(2) 



Pb 



(27) 



which is the familiar HNC approximation for the density 
profile p(r) in terms of the external potential and the 

(2) 

bulk direct correlation function c b (r) = c(r). Given c(r) 
from the previous HNC calculations of the pair structure 
in the bulk, Eq. ( p7| ) may be solved iteratively for any 
4>(r). If c(r) is replaced by its RPA form (^), then Eq. 
( p7| ) reduces to the MF form: 



p(r) = p b exp 



-(3<j>{r) 



Pb J dr>(|r-r'|) ( 



P(r') 
Pb 



- 1 



(28) 



Eq. (|2S]) also follows directly from the standard mean 
field approximation(MF-DFT) for for the intrinsic free 
energy functional |lq |: 

f3F n [p{v)]=f3F d + l - J drdr>(r,r>(r)p(r'), (29) 

which, in a different context, is identical to the functional 
used to derive the Poisson-Boltzmann theory for ionic flu- 
ids if v(r, r') is taken to be the Coulomb potential and p 
the charge density. 

Specializing to the case of a planar wall coinciding with 
the x — y plane, and confining the particles to the z > 
half-space, without any additional external potential, we 
note that the density profile p(z) satisfies the contact 
condition (l9|: 



p(z=0) = PP 



(30) 



where P is the pressure exerted by the particles on the 
wall, equal to the bulk pressure in the absence of an 
external potential. The sum rule |I(]) is satisfied by 



the MF-DFT approach, where pMF(0)=f3PMF=PbZMF, 
with Zmf defined by Eq. (0) since Zmf=Zh PA . How- 
ever, the sum rule is not satisfied by the (more accurate) 
HNC approximation (^?j) , which instead leads to (20) : 



P(0) = 2» 



\9pb ) ■ 



(31) 
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FIG. 8. Ratio phnc(0)/p(0) from ( |3l| ) compared to 
the analytic ratio obtained from the RPA2 approxima- 
tion to Eq. ([H])) for e* = 2. Also included is the ratio 
Pmf(0)/p(0) ~ Z C RPA /Z V HNC . Clearly the HNC is the 
better approximation, even though it does not satisfy the 
sum- rule of Eq. @ 



This reduces to the exact result ( |30[ ) provided the pres- 
sure is a quadratic function of the bulk density. This is 
very nearly true over a wide range of densities, as shown 
in the previous section. In particular the simple RPA 
compressibility (or MF) e.o.s. (p^|), which provides a fan- 
representation of the numerical HNC results, is of the 
necessary linear form to make Eq. ( p?i| ) and Eq. ( |3l| ) 
compatible. The deviations from the sum-rule (30) may 
be traced back to the slight non-linearity of ZhnC: as 
demonstrated in the inset of Fig. || and in Fig. 0, The 
relative error does not exceed 3% for e* — 2, although 
it tends to increase with increasing e* . In fact, the ratio 
Phnc{Q)I p(0) may be estimated from the very accurate 
RPA2; since Z P p A 2 — Z^, P A , the required pressure may 
be calculated from Eq. ( |23|) , while the RPA2 inverse com- 
pressibility is calculated to be : 



fd[3P\ 



2a 



Li 1/2 (-a) - Li 3/2 (-aj] 



(32) 



The resulting analytic estimate of Phnc(0)/p{0) is also 
shown in Fig. p|; as expected, it gives a good approxima- 
tion of phnc(0)/p(0). Even though the MF approach 



G 



exactly satisfies the sum-rule (pO|), the HNC approach, 
which does not satisfy (|3C|), is a better approximation. 
We note that the arguments above can be extended 
to the popular Percus-Yevick approximation |jT3| , where 
p PY (Q)/p(0) ~ (l + l/2a)<-l/2) Q, which, in contrast 
to the HNC or MF approaches, becomes increasingly less 
accurate as the density increases. 

We have numerically solved the HNC and MF Euler- 
Lagrange equations (|27j) and (|2^) to calculate the den- 
sity profiles p(z) of a GCM near a hard wall, for sev- 
eral values of the bulk density The theoretical profiles 
are compared in Fig. [)]to the results of MC simulations. 
The agreement is seen to be excellent, particularly at the 
higher densities. In fact, within the accuracy of the fig- 
ure the difference between the HNC and MF approaches 
is visible only for small z at p = 0.1, where, as expected, 
the HNC approach is slightly more accurate. 



P(z)/P 




FIG. 9. Density profiles from HNC and the MF-DFT 
for Gaussian particles (e = 2, R = 1) near a hard wall. 
Symbols are for MC simulations at 3 densities, the solid 
lines are from the MF-DFT approach, and the dashed 
lines are from the HNC approach. The two theories and 
simulation agree to within the accuracy of the graph for 
p* = 0.5 and p* = 1, but small discrepancies appear for 
p* = 0.1, where the HNC is slightly more accurate. 



In Fig. [K] we show density profiles for particles inter- 
acting with an external potential f3<fi(z) = exp[— z]/z, a 
situation similar to that encountered for polymer coils 
near a wall ||. Once again, we observe that the HNC 
and MF-DFT approaches are very close, implying that 
both are very accurate and could potentially be fruitfully 
combined with the effective potentials between polymer 
CMs H to derive a full DFT for polymer solutions in 
complex geometries. 

In summary then, the results of this section confirm 
that the model considered indeed behaves as a "mean 
field fluid" under inhomogeneous conditions. 
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FIG. 10. Density profiles from HNC (symbols) and 
MF-DFT (lines) for Gaussian particles (e = 2, R = 1) in- 
teracting with an external potential fl4>(z) = exp[— z]/z. 



IV. PHASE SEPARATION IN 
TWO-COMPONENT REPULSIVE GAUSSIAN 
MIXTURES 

Since the underlying polymer mixtures exhibit inter- 
esting phase behavior under a variety of physical condi- 
tions, it is natural to consider binary mixtures of Gaus- 
sian core mixtures, interacting via pair potentials: 



Pv Vil {r) = e* expHr/i^) 2 ], 



(33) 



where the species indices 1 < v, p, < 2. The total num- 
ber density is still denoted by p = (Ni + N2)/V , while 
the concentration variable x — N2/N. We are inter- 
ested in the possibility of a phase separation, or de- 
mixing transition, of the two species. The thermody- 
namic stability conditions for any binary mixture can be 
expressed in terms of the Helmholtz free energy per par- 
ticle, f(x, v) — F(Ni,N 2 , V)/N, considered as a function 
of the intensive variables x and v (or p — 1/v), for any 
fixed temperature. These conditions are [Bjl: 



> 



dv 2 



dx 2 



dv 2 
( d 2 f \ 2 



> 0. 



(34a) 
(34b) 
(34c) 



The first inequality expresses mechanical stability, (i.e. 
positive compressibility), the second is the condition for 
stability against spontaneous de-mixing at constant vol- 
ume while the last inequality ensures stability at constant 
pressure; it is equivalent to the more familiar condition 
{d 2 g(x,P)/dx 2 )p > 0, where g is the Gibbs free energy 
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per particle. We note that if ei ther of the first two sta- 
bility conditions (34a) or ( 341;) ar e violated, the more 
restrictive stability condition, (34c), is v iolated as well. 
Spinodal instability occurs when Eq. ( |34c| ) , is satisfied as 
an equality. The condition is equivalent to the k — ► 
divergence of the concentration-concentration structure 
factor: 

S cc {k) = x 2 S 11 (k) + (1 - x) 2 S 22 (k) - 2x(l - x)S 12 (k), 

(35) 

where the S Vfi (k) are the usual partial structure factors 
|jl3|| . From the OZ relations for a binary mixture it is 
easily inferred that S cc (0) diverges when: 



[l-(l-ar)pcii(0)] [l-xpc 22 {Q)\ 
-x(l-x)p 2 [c 12 (0)] 2 =0. 



(36) 



We now examine the implications of these conditions 
within the MF approximation, which we have shown to 
yield reliable results, except at low reduced density p*. 
The MF free energy ([l9]) , properly generalized to the bi- 
nary situation, reads: 

f(x, P ) = r d { P ) + r ix {x) + l - P %{x) (s?) 

where the first, ideal gas, term is irrelevant in the subse- 
quent considerations, f mlx is the ideal mixing term: 

f mix (x) =x\ilx + (1-x)Iil(1-x), (38) 

and the MF interaction term is: 

V (x) = (1 - x) 2 /% n (0) + 2x(l - x)[3v 12 {0) + x 2 /3v 22 (0). 

(39) 



The {/30„ M (0)} are the k -> limits of the FTs of the 
interaction potentials. In fact, the MF free-energy j37| ) 
has the same mathematical form as a second- virial theory 
which would be valid for very low densities p2[ . 

W it h th e MF free energy (^7|), the stability conditions 
(34a - 34c) reduce to: 



l + pV (x) > 
1 — px(l — x)x > 
1 + pVl(x) -p 2 x{l-x)A >0 



(40a) 
(40b) 
(40c) 



respectively, where the following parameters were de- 
fined: 



X = 2(3v 12 (0) - (/?fi u (0) + /3t) 22 (0)) 
A = ((3v 12 (0)) 2 -p0 u (0)/3v 22 (0) 



(41) 
(42) 
(43) 



Eq. (40a) can only be violated if the potentials them- 
selves violate a 2-component extension of Eq. ( A5) from 



Appendix A, which is a necessary (but not sufficient) con- 
dition for the existence of a well-defined thermodynamic 
limit. The limit of stability of the mixture (i.e. the spin- 
odal line) at cons tant volu me or pressure is reached when 
the inequalities (40b) and (40c) turn into equalities; the 



latter condition also follows from (|36|), when the c„ M are 
replaced by their RPA limits c VfJ .(k) — —j3v v ^(k) . De- 
mixing at constant volume is possible, provided x > 0; 
the density along the spinodal is then easily calculated 
to be: 



Pa{x) 



1 



x(l - x)x 



(44) 



De-mixing at constant pressure is only possible provided 
A > 0. The corresponding density along the spinodal 
satisfies: 



Vl(x) + JVi(x) 2 + M1 - x)A 
P*( x ) = 3a > ( 45 ) 



2x(l - x)A 
the pressure along the spinodal is: 

P s (x) = p s {x) + -p 2 (x)V (x), 



(46) 



and the critical consolute point is determined by the con- 
dition: 



dP s (x) 
dx 



= 0. 



(47) 



The simple quadratic expression for the pressure P (|4q), 
is easily inverted to obtain an expression for the spinodal 
density as a function of concentration x and pressure P. 

We now apply these general considerations within the 
MF framework to the binary Gaussian core model for 
which 



= n 3/2 e* Ufl R 3 Ufl . 



(48) 



Inserting the binary GCM expression for {^(O) into ex- 
pressions d42| ) and ( p[3| ) for x and A, we find that phase- 
separation at constant volume or at constant pressure is 
possible provided: 

X = tt 3 / 2 [2e* 12 R\ 2 - (e* n R 3 n + e* 22 R 3 22 )] > 0. (49) 



A = n 3 [(e* 12 ) 2 Rf 2 - e* n e* 22 R 3 n R% 2 )] > 0. 



(50) 



In order to focus on physically relevant values of the 
parameters e* M and R vfl , it is important to make contact 
with known results for polymer coils in a good solvent 
^,0. Simulations on binary solutions of self-avoiding 
polymer coils carried out in the low concentration limit 
M suggest that the effective pair potentials between cen- 
ters of mass are reasonably well represented by the Gaus- 
sian form (E53|), with : 



-11 — e n 



(51) 



8 



and 



15 



R 



12 



R 



22) 



(52) 



The relation (51) between the e ufl favors mixing. On the 
other hand i?i2 > l/2(i?n + -R22), which resembles the 
positive non-additivity that can drive de-mixing in hard- 
core mixtures 24 . Substituting ( |52| ) into (|49|), we find 
that a spinodal instability of the mixture is possible at 
constant volume provided: 



> sfl- 



1 + (R22/R11) 



(1 + (iWi?n) 2 ) 



3/2 



> 1 



(53) 



10 - 



pRn 



FIG. 11. Constant pressure spinodal ( |45| ) for param- 
eters taken from simulations of L = 100 and L — 200 
monomer effective polymer CM potentials [^| . The x-axis 
denotes the composition x = X2 = N2/N. The y-axis de- 
notes the density pRfi, where Rn is the = radius of 
gyration Rg for the L — 200 polymers. The dot is the 
critical point at (x = 0.70, pRii = 5.6). 



which contradicts the requirement (|51|). On the other 
hand, if (|||) is substituted into (|50|), de- mixing at con- 
stant pressure may occur provided: 



> 



2{R 22 /Ru) 



1 3 



1 + (R22/R11) 2 



< 1, 



(54) 



which is compatible with the requirement (pr]). 

More specifically, we have chosen values of the pa- 
rameters e* M and i?„ M appropriate for a polymer mix- 
ture of self- avoiding polymers of L = 200 (species 1) and 
L = 100 (species 2) monomers |f^,p3|. The resulting spin- 
odal line ( |45| ) in the x — p plane, calculated from the MF 
free energy (|37|) with (|48|), is shown in Fig. 
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Phase 

separation into two solutions of different composition x 
occurs above a critical density p* = 5.6/ R^ and critical 
composition x c — 0.70. Note that since all terms in the 
free energy are of entropic origin, the temperature scales 
out, i.e. the mixture behaves as an athermal system. In 
view of the remarkable accuracy of the MF theory at high 
density, as illustrated in sections II and III for the one 
component GCM fluid, we expect the phase-diagrams, 
calculated within MF (or equivalently RPA) to be reli- 
able; full calculations of the binodal line, based on RPA2 
and HNC theories, will be reported elsewhere. 



V. CONCLUSIONS 

The calculations carried out in this paper, and in re- 
lated work lead to the conclusion that a system of 
classical particles interacting via a repulsive Gaussian 
core potential behaves like a weakly correlated "mean 
field fluid" over a wide range of temperatures and densi- 
ties. In fact for any temperature there is always a (sur- 
prisingly low) density beyond which the excess free en- 
ergy per particle is a linear function of density and the 
resulting excess pressure increases like a quadratic func- 
tion of the density. On the other hand, in the opposite 
low density regime, a virial expansion of the equation of 
state in powers of the density appears to converge only at 
extremely low densities. This is in sharp contrast to hard- 
core systems, for which the virial expansion provides a 
good estimate of the equation of state up to relatively 
high packing fractions |17|| , while the pressure diverges 
near close packing according to a simple free volume pic- 
ture. At very strong interaction strength (e* > 100), 
the GCM behaves effectively as a hard core fluid that 
freezes at intermediate densities, but re-melts under fur- 
ther compression to return to mean-field like behavior 
P,^0[. The small correlational effects at low and inter- 
mediate densities are adequately described by the simple, 
analytic RPA2 extension of RPA theory, or by the HNC 
integral equation (requiring numerical solution) , which is 
nearly thermodynamically consistent over a broad range 
of temperatures and densities. 

The MF theory performs equally well in the inhomo- 
geneous situation of Gaussian core particles near a hard 
wall. The binary version of the model phase-separates at 
high densities, when the widths of the Gaussian repulsion 
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satisfy the composition rule (p2|), provided condition (KJ) 
is satisfied. This provides an interesting example of phase 
separation in systems with purely repulsive interactions. 

To conclude it seems worthwhile to consider the rel- 
evance of the GCM for the description of polymer so- 
lutions. The latter enter the semi-dilute regime when 
polymer coils start to overlap, i.e. when p* ~ 3/(4tt). 
For densities of this order we have seen that the GCM 
behaves like a "mean field fluid" , with a quadratic density 
dependence of the pressure. The exponent 2 is close to 
the 9/4 power observed for the osmotic pressure of semi- 
dilute polymer solutions [£5| . The difference between the 
exponent 9/4 and 2 is due in part to the weak, but sig- 
nificant density dependence of the effective pair potential 
between the centers of mass of self-avoiding polymers || , 
which leads to an additional density dependence of the 
RPA or MF equation of state (O). This possibility is 
being explored in more detail |26|. 

The effective polymer- wall potentials derived in rcf. ||] 
show a significant variation with density. Nevertheless, 
the form of the p(z)/pb for the GCM in a fixed external 
potential follows the same qualitative trends as the dis- 
tribution of the polymer CM's near a wall, pCM(z)/pb, 
suggesting that the physics of polymer coils near a wall 
is well captured by the GCM. 

The de-mixing transition of binary Gaussian core mix- 
tures is reminiscent of the tendency of polymers of differ- 
ent molecular weight to phase separate at high concentra- 
tion and in the melt. Again, further analysis is required 
to decide if the analogy between the de-mixing of Gaus- 
sian core mixtures and of polymer blends is fortuitous, 
or has some deeper foundation. 
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APPENDIX A: THERMODYNAMIC STABILITY 
OF SOFT-CORE POTENTIAL SYSTEMS 



It was pointed out in section II A that the GCM sat- 
isfies Ruelle's condition (||) for the existence of a finite 
thermodynamic limit. In this appendix we give two ex- 
amples of pair potentials involving a repulsive core and 
a small attractive component which do not satisfy Ru- 
elle's stability condition, and hence belong to the class of 
potentials referred to by him as "catastrophic" . The fol- 
lowing considerations are not completely academic, since 



it has been shown in ref |2(| that the effective pair po- 
tential between the centers of mass of two polymer coils 
in a good solvent indeed exhibits a small attractive part 
at distances of the order of several times the radius of 
gyration R g for intermediate densities. When the poly- 
mer coils are no longer in a good solvent the potentials 
can develop even larger attractive parts . 

According to proposition 3.2.2 in Ruelle jl^], given an 

(2) 

interaction energy built up by pair potentials vi , 
the grand partition function is finite only if the follow- 
ing two equivalent properties hold for all N > and all 
{r ; } e R N : 



N N 



(Al) 



and 

4 2) (r 1 ,....,r N )= Yl v 2 (\r l -r 3 \)>-NB (A2) 



l<i<j<N 



for a B > 0. Note that in (Al) the double sum includes 
the self-interaction (i — j). 




□ p=0 polymer potential 

V A (r): catastrophic potential 

V Ei( r ) : catastrophic potential 



FIG. 12. Two "catastrophic" potentials compared to 
a typical CM potential for two polymers in a good sol- 
vent. _Potential va{t) (A3) violates the conditions ( |A l| ) 



and (A2) for h omo geneous fluid configurations, while po- 



tential vb(t) (A4) violates (Al) and (A2) only for inho- 



mogeneous configurations like the fee crystal. 

We consider two examples of potentials which do not 
satisfy these conditions: 



va{t) — 1.87 cos 
vb{t) = 1.87 cos 



V(2 + £) 



1.7R G 



1.7R G 

exp 



exp 



1.7 Rg 

4' 



1.7R G 



and compare them in Fig. [IJ to the polymer CM po- 
tential between two isolated L = 500 SAW polymer 
coils. Here S is an arbitrary positive constant taken to be 
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8 = 0.001 in Fig. |g Although at first sight they don't 
appear very different from the purely repulsive polymer 
potential, they are both "catastrophic" . 

If S > 0, the first potential, Va{t), violates a weaker 
condition than (Al) or (A2), namely 



«(0) 



v{r)dr > 0, 



(A5) 



which is necessary (but not sufficient) for a thermody- 
nami c li mit. When it doesn't hold, conditions (Al) 
and (A2) can be violated for a homogeneous "gas" with 
g(r) — 1. This has a further implication for fluids de- 
scribed by a mean-field free-energy ( |l9| ) , since the inverse 
compressibility dflP/dp = 1 + f3v(0)p cannot go through 
zero without violating the condition flA5|), which implies 
that one-component soft-core fluids described by a mean- 
field e.o.s. cannot support a spinodal instability. 

The second poten tial, Vb(t) , ha s an integral £>_b(0) > 0, 
but it still violates (Al) and (A2) for an inhomogeneous 
configuration. For example, for an fee lattice with single 
occupancy £\ £\ v B (\rj - r ; |) = -0.13iV. The po- 
tential is catastrophic because one can always lower the 
total energy indefinitely through multiple occupancy of 
the lattice sites. 

The 9 point in polymer solutions can be defined as the 
temperature where the effective second osmotic virial co- 
efficient, £?2, passes through 0. Above the 9 point the 
solvent is said to be "good" , while below the 9 point the 
solvent is said to be "poor" . Simulations of a model for 
two polymers in a poor solvent show that the effective 
pair potential is no longer strictly positive definite be- 
low the 9 point |tJ , implying that the pair potentials can 
become catastrophic. In fact, for the type of polymer 
CM potentials considered, this seems to occur just below 
the 9 point temperature where B 2 = 0. It is tempt- 
ing to speculate that the coil-globule transition, which 
also typically occurs slightly below the 9 temperature, 
is related to the point at which the effective pair poten- 
tial becomes catastrophic. However, it is not yet clear 
whether the pair-potential picture of polymer solutions 
H remains valid for poor solvents. 



APPENDIX B: VIRIAL EXPANSION FOR THE 
GCM FLUID 

In this appendix we briefly consider the convergence 
of the virial expansion of the equation of state of the 
GCM in powers of the density p. The FT of the Mayer 
f-function in Eq. ( [Toj ) is given by the convergent sum: 



/(*) 



3/2' 



ex P(-3^)(- £ 



n=l 



(Bl) 



Here the width parameter R in the Gaussian potential 
(|l|) has been chosen as unit of length for convenience. 
The second and third virial coefficients B 2 and B3, of 



the GCM can then be expressed as the following conver- 
gent sums: 



Bo 



1 Tr 3 / 2 °° 

>)-VE 



(-£*)" 



2 nln 3 / 2 

n—l 



00 00 00 



I ~- ' ~- I c *\i+j+k 

B = 1 3y^x^Y^ l~ £ ) 

3 3 l^l^l^ i\j\k\(ij + jk + ikW 2 
i=l J=l k=l J ' 



(B2) 
(B3) 




FIG. 13. Second virial coefficients for a Gaus- 
sian potential as a function of interaction strength e*. 
(-£?! > B2 > E>2 lr ). Also included is the empirical rela- 
tion B| /S (v' ln (2 e *)) 5 where B| /S (cr) is the hard-sphere 
second- virial coefficient. 

The variations of B2 and B3 with e* are shown in Figs 
|l3| and |l4|; both virial coefficients are always positive. 
The virial expansion of the equation of state reads: 



(3_P 
P 



Z = — = 1 + B 2 p + B 3 p 2 + 0{p 3 ) 



(B4) 



and the results from the 2 and 3 term series are compared 
in Figs. |^-^ of section |H| to the predictions of the RPA 
and HNC theories, and to MC simulations for e* = 2, 10 
and 90. The virial expansion is seen to break down very 
early. In particular, although the MF e.o.s., which be- 
comes very accurate at high density, predicts a linear 
variation of Z with density, the slope differs more and 
more from B2 as the interaction strength e* increases. 
Adding the B 3 contribution leads to rapid deterioration 
of the predicted e.o.s. as the density increases. 

The shortcoming of the virial expansion in powers of 
density is further illustrated by considering the RPA. 
From Eqns. ( jig ) and (^3|), one may extract the following 
compressibility and virial estimates of the 2nd and 3rd 
virial coefficients: 



B^ = -^' 2 e*;Bl 
2 2 . 3 







gvir _ 1^3/2^* J ^ _ I . nvir 



7T 3 (e*) 3 
9\/3 



>^3 



(B5) 
(B6) 
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As shown in Figs. [13] and [LJ, the exact virial coefficients 
are bracketed by the virial and compressibility estimates 
extracted from the RPA: 



for n > 3. This, B^ r (B6) is added, recovers exactly 
the virial RPA equation of state (^3j), which can only be 
expanded in powers of density for a < 1, i.e. for 



B c 2 > B 2 > B% lr 
BY >B 3 >B C 3 = Q 



(B7) 
(B8) 



The large deviations shown in Figs. 13 and 14 imply that, 



in contrast to the case at high densities, the RPA is ex- 
pected to perform poorly at very low densities and large 
e*, where it is thermodynamically inconsistent. In fact, 
BV? r even goes negative for e* > 5.7! The effect this has 
on the RPA virial e.o.s. is demonstrated in the inset of 
Fig. H and in Fig. |l5|. However, even though for e* = 2, 
Bf r is 13% less than B 2 , and Bf r is over 300% larger 
than £?3, Z V RPA remains within 1% of the exact e.o.s. over 
the entire density range! Thus, in spite of the fact that 
the RPA virial approximation grossly misrepresents the 
first two virial coefficients, it nevertheless accurately de- 
scribes the e.o.s., implying that the density is not a good 
expansion parameter for the GCM fluid phase. 
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FIG. 14. Third virial coefficients for a Gaus- 
sian potential as a function of interaction strength 
e*. B% = 0, Bf r > B 3 . Also include d is the 
empirical relation Bs(e*) 



B? s 



Bf s (V ln(2e*)), where 
(a) is the hard-sphere third-virial coefficient. In- 

3/2, 



set: p m = 1/(tt j/ V) is the maximum density for which 



the RPA virial e.o.s. ( p 

sion in powers of the density. p cp — \/2/(a^/ f ) 3 is the 



can be written as an expan- 
Pcp = V2/(<r?f S f) 
density at which the effective hard-sphere system with 
the same second virial coefficient as the GCM would be 
close-packed, and beyond which a virial equation would 
no longer be expected to exist. 

A further hint at the breakdown of the virial expan- 
sion comes from summing the virial series to all orders in 
the high temperature limit, where, from the diagramat- 
tic representation of the virial coefficients [|l3| , it can be 
shown that the B n are given by: 



1 ac* 

Bn = --7T 



I) (-e*)"(n- 1) 



7,5/2 



0(e*) 



n+l 



(B9) 



P < Pn 



r 3/2. 



0.1796/e* 



(B10) 



This implies that in the high-temperature limit, the virial 
expansion does not converge for densities higher that p* m . 
For e* = 2, there is no convergent density expansion of 
the RPA virial e.o.s. for p* > p* n w 0.0898. A similar 
breakdown in convergence may be expected for the ex- 
act virial expansion. The physical reason for this lack 
of convergence lies in the possibility of multiple overlap 
of soft core particles, giving much more weight to higher 
order cluster integrals compared to the case of fluids with 
hard-core interactions. 

At large enough e* , the overlap probability becomes 
exponentially small, and the GCM can be mapped onto 
an effective hard-sphere system ||[l(|. One possible cri- 
terion for the mapping is to equate the second virial co- 
efficients. From this we obtain an effective hard-sphere 
radius of: 



HS 
T eff 



2tt 



B 2 (e* 



1/3 



(Bll) 



which for e* > 1 is well app roximated by the empirical 
expression cr^A ~ ^/ln(2e*). For large e* and low densi- 
ties, the equation of state resembles that of hard-spheres 
(see e.g. Figjl5|), suggesting that a virial expansion does 
indeed exist for low densities. We note that for this large 
value of e*, the true virial expansion appears to have a 
larger radius of convergence than that of the RPA virial 
e.o.s., for which p m {o^ f s f f « 0.023. For e* > 100 there is 
a freezing transition at roughly the density expected for 
the effective hard-sphere system (p(crfy-y ) 3 w 1), not far 
above which any effective virial expansion is expected to 
break down (see e.g. the inset of Fig. Q). In fact, since 
Gaussian potentials don't have an infinitely hard core, 
it is possible to achieve much higher densities than are 
normally available to simple liquids. At the lowest densi- 
ties the fluid is described by a linear second virial theory 
e.o.s., but as the density increases, this rapidly turns over 
to a mean field like linear e.o.s. with a different (larger) 
slope. Thus, even though the e.o.s. is well described by 
a first order polynomial in the density p it is not at all 
equivalent to a second virial theory, and the density is 
generally not a good expansion parameter. 
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FIG. 
Here e* 



15. Z v.s. p(a e ff) in the low density limit. 
= 90 so that (a^/f) = 2.27; p(°" S f ) 3 = 1 corre- 
sponds to p* — 0.085. For low effective density the e.o.s. 
follows the hard-sphere e.o.s. (here approximated by the 
Carnahan-Starling form [Q). For higher densities the 
fluid moves towards the mean-field fluid limit (see Fig. 
Q). Note that the two RPA expressions for the e.o.s. are 
very poor approximations in this low density regime. 
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